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Abstract 



Pricing exotic multi-asset path-dependent options requires extensive Monte Carlo sim- 
ulations. In the recent years the interest to the Quasi- monte Carlo technique has been 
renewed and several results have been proposed in order to improve its efficiency with 
the notion of effective dimension. To this aim, Imai and Tan introduced a general 
variance reduction technique in order to minimize the nominal dimension of the Monte 
Carlo method. Taking into account these advantages, we investigate this approach in 
detail in order to make it faster from the computational point of view. Indeed, we real- 
ize the linear transformation decomposition relying on a fast ad hoc QR decomposition 
that considerably reduces the computational burden. This setting makes the linear 
transformation method even more convenient from the computational point of view. 
We implement a high-dimensional (2500) Quasi-Monte Carlo simulation combined with 
the linear transformation in order to price Asian basket options with same set of param- 
eters published by Imai and Tan. For the simulation of the high-dimensional random 
sample, we use a 50-dimensional scrambled Sobol sequence for the first 50 components, 
determined by the linear transformation method, and pad the remaining ones out by 
the Latin Hypercube Sampling. The aim of this numerical setting is to investigate the 
accuracy of the estimation by giving a higher convergence rate only to those components 
selected by the linear transformation technique. We launch our simulation experiment 
also using the standard Cholesky and the principal component decomposition methods 
with pseudo-random and Latin Hypercube sampling generators. Finally, we compare 
our results and computational times, with those presented in Imai and Tan [8]. 

Key Words: Effective dimensions. Path-generation techniques. Linear transforma- 
tions. Quasi-Monte Carlo simulations. 

1 Introduction 

The Monte Carlo method (MC) is a computational intensive technique whose purpose 
is to estimate integrals numerically. It is characterized by a rate of convergence of order 
0{1/ -Jn), where n is the number of simulations, and it is independent of the problem 
dimension d. This last feature makes the MC method appealing and applicable to 
several financial high-dimensional situations such as options pricing. Furthermore, the 
estimation error (RMSE), that can be easily computed statistically, depends only on 
the convergence rate and on an intrinsic constant. 

Based on probabilistic considerations, standard reduction techniques can only reduce 
the constant but cannot improve the convergence rate. 

In contrast, Quasi-Monte Carlo methods (QMC) aim to enhance the convergence 
rate by means of low-discrepancy sequences. These sequences provide better stratifi- 
cation and a convergence rate of order O {^^77^ ) (see Niederreiter [ID]). The rate is 




faster than the previous one but depends on the problem dimensions. These sequences 
are purely deterministic, meaning that the estimation error cannot be estimated statis- 
tically. In the Randomized Quasi-Monte Carlo (RQMC) method some randomness is 
introduced in the low-discrepancy sequences while preserving their better convergence 
rate. This technique is called scrambling. 

Several numerical investigations conclude that QMC and RQMC simulations do not 
give substantial advantage for d > 10/20. 

Some approaches have been proposed in order to extend the QMC superiority to 
high-dimensional estimations. Caflisch et al [T] address the problem using the analysis 
of variance (ANOVA) of the integrand function and defining two notions of effective 
dimension: the effective dimension in truncation and superposition sense. Briefly, the 
truncation dimension reflects that, for some integrand functions, only a small number of 
inputs really matter. The definition of effective dimension in superposition sense takes 
into account that for some integrands the inputs might influence the outcome through 
their joint action within small groups. 

Imai and Tan [B] proposed a general linear transformation construction (LT) to 
reduce the effective dimension of the problem in superposition sense, focusing on the 
particular payoff function. The authors show that this approach offers a considerable 
advantage with respect to the principal component analysis (PCA) in terms of accuracy 
and versatility. 

Moreover, their simulation procedure relies on the complete Latin Supercube Sam- 
pling generation (see Owen [llj for more on this topic) in order to generate a high- 
dimensional low discrepancy sequence with good properties. 

Here we investigate the accuracy of the LT method in detail and implement the 
construction fast by an efficient QR decomposition. We run our simulation procedure 
with the same set of parameters as in Imai and Tan fSl and are thus able to directly 
compare the respective results. 

We will demonstrate, that our implementation makes the LT considerably faster 
and maintain its versatility. 

We test the efficiency of the LT construction by launching a MC simulation in 
a more extreme setting. We use scrambled low-discrepancy sequencies only to those 
components the LT considers as optimal, while simulating the others with the Latin 
Hypercube Sampling (LHS) that has lower convergence rate. 

The LHS is supposed to give good accuracy when the target function is a sum of 
one-dimensional ones. If the LT accomplishes this task optimally it would give a good 
improvement in this setting too. Our experiment is intended to test if the LT gives the 
same results as in Imai and Tan [H], in terms of RMSEs, in this partial RQMC setting. 
This means, that if the LT with RQMC provides considerable advantage with respect 
to the pure LHS generation it reduces the effective dimension in superposition sense 
really optimally. 

As a comparison, we launch the MC simulation using a standard pseudo-random 
generator and build the random path with standard Cholesky and PCA decompositions 
too. 

The paper is organized as follows. Section 2 describes the financial setting and 
formulates the Asian basket option pricing problem as an integral explicitly. Section 
3 introduces the MC and the QMC methods and the notion of effective dimensions of 
the problem. Section 4 describes the LT construction introduced by Imai and Tan and 
how it applies to several financial situations. Section 5 presents the main steps of our 
MC simulation. Section 6 illustrates the numerical results we obtain and discuss the 
efficiency of the LT and its fast implementation. Section 7 concludes the paper and the 
Appendix describes the ad hoc QR decomposition used. 

2 Problem Statement 

We consider the problem of estimating the fair price of a contract in a standard financial 
market 371 in a Black-Scholes framework, with a constant risk-free rate r and time- 
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dependent volatilities. There are M + 1 assets in the market, one risk free asset and 
M risky assets. The price processes of the assets in this market are driven by a set of 
stochastic differential equations. 

Suppose we have already applied the Girsanov theorem and found the (unique) risk- 
neutral probability, the model for the risky assets is the so called multi-dimensional 
geometric brownian motion: 

5oW = e^* (1) 
dS,{t)^rS,{t)dt + a,{t)S,it)dW,it), i = l,...,M. (2) 

Here Si {t) denotes the i-th asset price at time t, ai (t) represents the instantaneous time- 
dependent volatility of the i-th asset return, r is the continuously compounded risk-free 
interest rate, and W (t) = {Wi {t) , . . . , Wm {t)) is an M-dimensional Brownian motion. 
Time t can vary in M^, that is, we can consider any maturity T G R*^ for all financial 
contracts. 

The multi-dimensional brownian motion W [t) is a martingale, each component is 
a martingale, and satisfies the following properties: 

E[Wi{t)]^0, i = l,...,M. 

[W^,Wk]{t)^ p,kt, i,k^l,...,M. 

where [vK^) represents the quadratic variation up to time t and pik the constant 
instantaneous correlation between Wi and Wk ■ 

Applying the risk-neutral pricing formula, the value at time t of any European 
T-maturing derivative contract is: 

V{t)^exp{r{T^t))E[4>{T)\Tt]. (3) 

E denotes the expectation under the risk neutral probability measure and 4>{T) is a 
generic Tt measurable function, with Tt = cr{Q < t < T; W(t)}, that determines the 
payoff of the contract. Although not explicitly written, the function 0(r) depends on 
the entire multi-dimensional brownian path up to time T. 

We will restrict our analysis to Asian options that are exotic derivative contracts 
that can be written both on a single security and on a basket of underlying securities. 
Hereafter we will consider European-style Asian options whose underlying securities 
coincide with the M + 1 assets on the market. This is the most general case we can 
tackle in the market 971 because it is complete in the sense that we can hedge any 
financial instrument by finding a portfolio that is a combination of this M + 1 assets. 



2.1 Asian Options Payoff 

The theoretical price for a discretely monitored Asian option is: 



Qi (t) = exp {r{T - t)) E 



a{t) = exp{r{T -t))E 



N 



K 



M N 



Option on a Single Asset 

(4) 

Option on a Basket (5) 



where t\ < t2 ■ ■ ■ < tj^ = T and the coefficients Wij satisfy 

European options with payoff functions ^ and (O are called arithmetic weighted 
average options or simply arithmetic Asian options. When M > and iV = 1 the 
payoff only depends on the terminal price of the basket of M underlying assets and the 
option is known as basket option. 
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2.2 Problem Formulation as an Integral 



The model 9Jt, presented in the first section, consists of the risk-free money market 
account and M assets driven by M geometric brownian motion described by equation 
^ whose solution is: 



Si (t) = Si (0) exp 



ds 



a, (s) dW, is] 



,i = l,...,M. 



(6) 



t ^ f 1 

The quantity ds is the total volatility for the i-th asset. The solution ([H]) is a 



multi-dimensional geometric brownian motion, written GBM^r, ^-^^^ds 

Under the assumption of constant volatility the solution is still a multi-dimensional 
geometric brownian motion with the following form: 



S, (t) = S, (0) exp 



,i = l,...,M. 



(7) 



In compacted notation the solution ([7]) is GBM^^r, ^tj . 

Pricing Asian option requires to monitor the solutions © and ([7]) at a finite set of 
points in time {ti, . . . , tjf}. This sampling procedure yields to the following expressions 
for time-depending and constant volatilities: 



Si{tj) = St{0)exp 



Si{tj) = Si{0)exp 



dt + Z,{tj) 



(8) 



(9) 



where the components of the vector (Zi(ti), . . . Zi(tAr), ^2(^1), . . . , ZM{tN)) are M x N 
normal random variables with zero mean vector and the following covariance matrix: 



Smjv — 



for time-dependent volatilities, (10) 



for constant volatilities. 



Each element depends on four indexes: 



CTi {t)(7k {t)pikdt 



with i,k = 1, . . . , Af and l,m ~ I, . . . , N. 

The payoff at maturity T of the arithmetic average Asian option is then: 



where 



and 



Pa{T)^{g{Z)-K)+ 

MxN 

g{Z) = ^ exp {fik + Zk) 



k=l 



^fe = \\i{wk^k^Sk^{Q)) + r 
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(11) 



(12) 

(13) 
(14) 

(15) 
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for constant volatilities or 

f*"^ al (t)dt 

fik^ln{wk,k,SkAO)) + rtk,-^ (16) 

for time-dependent volatilities. The indexes ki and k2 are respectively ki = {k — 
l)fnodM,k2 = [{k — l)/M] + 1, where mod denotes the modulus and [.] the greatest 
integer less than or equal to x. 

The calculation of the price a{t) in equation ([5]) can be formulated as an integral on 
[0, 1]^^^ in the following way (see Dahl and Benth 2J and 3 ): 

a{t) = exp{r{T~t)) I {g{vi) - K)+ F^\vi)dM (17) 

J [0,1]^" 



3 Problem Dimension 

The main purpose of the standard MC method is to numerically estimate the following 
integral: 



/ /(x)dx. (18) 

I can be seen as E [/(f/)], the expected value of a function /(.) of the random vector 
U that is uniformly distributed in hypercube [0, 1]''. 

MC methods simply estimate / by drawing a sample of n independent replicates 
Ui . . . ,Un of U and then computing the arithmetic average: 

_ _ 1 " 

/ = /n = -^/(f/,). (19) 

The Law of Large Numbers ensures that /„ converges to / in probability almost 
surely and the Central Limit Theorem states that / — /„ converges in distribution to a 

normal with mean and standard deviation o j \/n with a = \J (/(x) — I)^ dx. The 
convergence rate is than 0{l/n) for all dimensions d. The parameter a is generally 
unknown in a setting in which / is unknown, but it can be estimated using the sampled 
standard deviation or root mean square error (RMSE): 



RAISE : 



\ 



1 " 2 

— E(/(t^»)-^«) ■ (20) 
i=i 



When the nominal dimension d of the problem of estimating the integral (|18p is one, 
there are standard numerical techniques that give a good accuracy when / is smooth. 
Considerable problems arise when d is high. 

We aim to estimate the fair value of the Asian option of equation (|17p with an 
high-dimensional Quasi MC simulation as formulated in p^ . 

QMC method relies on the construction of deterministic sequences, also known 
as low-discrepancy sequences, that cover the hypercube [0, 1)*^ uniformly. We define 
the quantity £>* — D* (Pi, . . . ,P„) as the star discrepancy. It is a measure of the 
uniformity of the sequence {Pn}„gN* ^ [O; l)'^ ^^'^ it must be stressed that it is an 
analytical quantity and not a statistical one. 

A sequence {Pn}„gN* called low-discrepancy sequence if: 

D:iP,,...,P..)^O^^y (21) 

The following inequality, attributed to Koksma and Hlawka, provides an upper 
bound to the estimation error of the unknown integral with the QMC method in terms 
of the star discrepancy: 

\I-I\<D:VHK{f). (22) 
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Vhk if) is the variation in the sense of Hardy and Krause. Consequently, if / has a 
finite variation and n is large enough, the QMC approach gives an error smaller than 
the error obtained by the crude MC method for low dimensions d. 

It is well known that QMC methods loose the better accuracy in high dimension. It 
is then fundamental to capture the most important (in statistical sense) components or 
to reduce the nominal dimension of the problem by means of ANOVA considerations. 

Let A = {1, . . . ,d} denote the set of the independent variables for / on [0, 1]''. / 
could be written into the sum of orthogonal functions each of them defined in a different 
subset of A, that is depending only on the variables in each of these subsets: 

/(x) = /"W (23) 

uCA 

Now let |u| denote the cardinality of u and cr^ ~ J(/(x) — dx, cr^ — J /jj(x)^dx, 
(jg — 0, supposing a < +oo and |u| > it holds: 

= E ^« (24) 

uCA 

Equation (|24p partitions the total variance into parts corresponding to each subset 
u C A. The /„ enjoys some nice properties: if j G u the line integral Jjp ^ /u(x) dxj = 
for any Xk with k ^ j, and if m 7^ w / fu{x)fv{x) dx = 0. 

Exploiting the ANOVA decomposition, the definition of effective dimension can be 
given in the following ways: 

Definition 1. The effective dimension of f, in the superposition sense, is the smallest 
integer ds such that J2o<\u\<ds - P'^^ ■ 

The value ds depends on the order in which the input variables are indexed. 

Definition 2. The effective dimension of f , in the truncation sense, is the smallest 
integer dx such that J2uC{i dt} '^u ^ P'^^ ■ 

< p < 1 is an arbitrary level; the usual choice is p = 99%. 

The definition of effective dimension in truncation sense reflects that for some in- 
tegrands, only a small number of the inputs might really matter. The definition of 
effective dimension in superposition sense takes into account that for some integrands, 
the inputs might influence the outcome through their joint action within small groups. 
Direct computation leads to: ds < dr < d. 

4 Linear Transform Construction 

Imai and Tan ^ proposed a general LT method for path generation with main purpose 
to minimize the effective dimension in truncation sense of a simulation problem. 

The LT approach provides the same results as the PCA-based one, moreover proves 
to be more accurate and versatile in certain situations. 

Many studies demonstrate that the QMC pricing of certain speciflc derivative con- 
tracts is not substantially improved by the brownian bridge construction. This suggests 
to focus the attention onto the particular payoff function while even the PCA approach 
is applicable only for multi-dimensional normal random variables. In contrast, the LT 
generation focuses on the particular payoff function instead of the multi-dimensional 
brownian path. 

This method provides the best results for linear combinations of normal random 
variables. Imai and Tan [B], [7] and [5] investigated the practical improvement of the 
LT method by running very high-dimensional simulations for European options, bonds 
pricing in different dynamics (see the cited references for more on this topic). 

A n-dimensional random vector Y with covariance matrix Sj, can be characterized 
starting from a vector of independent standard normal variables e by the following 
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transformation: y — Ce, with CC^ — T,y. Imai and Tan consider the foUowing class of 
LT as solution of the previous general problem: 

(jLT ^ (jCh^ (25) 

where C'-^^ is the Cholesky matrix associated to the covariance matrix of the normal 
random vector to be generated and A is an orthogonal matrix, i.e. AA^ = I. 

The optimum C^^ is obtained by optimally choosing A so that the effective dimen- 
sion in the truncation sense of the problem of interest is minimized. 

Maximizing the explanatory variability of a normal vector with covariance matrix 
E consists in finding the optimum orthogonal matrix A* by iteratively solving the 
following optimization problem: 

MN 

max||C't^|P= max V C^'^A.k (26) 

p=i 

subject to II A.kll — 1 and A.k • A*; = for z = 1, . . . , A; — 1 and k < n (A*; indicates the 
columns that have been already calculated) 

represents the fc-th column vector and C^^ the fc-th row vector of C^^; the 
same notation holds for all the matrices. Imai and Tan [8] proves that this procedure 
achieves the same results, in terms of explained variability, of the PCA decomposition 
of Smat- Indeed: 

max||C'iT||2 = (C^^'A.k)^C^'^A.k = A^^EA.k (27) 

Hence the optimization problem is similar to seeking the k-th principal component. 

Finding the optimal matrix A is equivalent of finding the optimal QR transformation 
of C with CC^ = E where R = [C'-^^)'^ and Q = Aiw the sense described before. 

The PCA decomposition provides the best solution for normal random vectors with 
Q = and R = A^/^ with V and A the orthogonal matrix of the eigenvectors and A 
the diagonal matrix of all the eigenvalues in decreasing order respectively. 



4.1 Special Cases 

As for linear combinations of normal random variables the LT approach minimizes the 
effective dimension in truncation sense. It is established from standard statistics that 
a linear combination of normal random variables is still a normal random variable with 
mean and variance that depend on the linear combination. It is than trivial that an 
integral problem with the nominal dimension d that involves a linear combination of d 
normal random variables has an effective dimension in superposition sense equal to one. 
The LT procedure returns this results in truncation sense as an optimization procedure. 

Let /(z) be a linear combination of d normal random variables /(z) = X^iLi WiZi, 
with z ^ N{^; E) and constants m;, i = I, . . . , d. If C denotes the generic decomposed 
matrix of S then the above function can be expressed as: 

d 

/(e) = ^afeefc + M-w (28) 

k=l 

where ak = C.k ■ w and e is a c?-dimensional vector of standard and independent normal 
random variables. Furthermore the total variance of / is: 

d 

^' = E"fc (29) 

i—k 

The truncation dimension is the smahest integer dx that satisfies: 

dT 

Y.al>pa^ (30) 

i=k 
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As with the LT approach the optimal C is C^-^ = C^'^A that leads to: 

afe = A.k-B k = \...,d (31) 

where B = (C'^'')"^w. Consequently, minimizing the effective dimension in the trunca- 
tion sense is equivalent to maximizing the variance contribution due to the first com- 
ponent a\ and obtaining A.i. Iterating this procedure and imposing the orthogonality 
condition we get the optimal matrix A. It can be proven, see Imai Tan [B] or [5], that 
the optimal solution for k — \ la: 

A\ = ±^ (32) 

while for k = 1,. . . ,d the column vectors can be arbitrary but must satisfy the or- 
thogonality condition. Substituting this results into equation (|30p we are left with 
ai = ±||B|| and = for A: = 2, . . . , d. The original function / can be written as: 

/(e) = Ai-w±||B||ei (33) 

This is the best possible scenario for the dimension reduction. The LT approach reduces 
any nominal d-dimensional problem involving a linear combination of normal random 
variables into a one-dimensional problem in truncation sense. This means that the LT 
method rearranges the linear structure of the function for the best possible reduction. 
Let us now consider the following function: 



/(e) = expl n + y] atek - K (34) 




with /i, ak and K constant. 

,f{x)^ can be considered the payoff function of a geometric average Asian option 
with strike price K and: 

tJ^^YnLiYfj^iWi, logSi{0) + (r - ^tj^ , = C.k • w . (35) 

Such a derivative contract is not traded but nevertheless serves to understand the 
computational problem. Indeed, performing the logarithm log{f{e) — K), we obtain a 
new function which is a linear combination of normal variates. Applying the results of 
the LT method for the previous example we showed that the nominal dimension MN 
of the new problem shrinks to one. Again this is not surprising because we know that 
the product of log-normal variates is still a log-normal variate. 

These examples highlight the main differences between the PCA decomposition and 
the LT methods. The former returns the best decomposition of the covariance matrix of 
a normal random vector in terms of variability of each component. The latter reduces 
the effective dimension of the problem focusing on the particular payoff function. It 
provides the best solution for linear combinations of normal variates. 



4.2 General Case 

General payoff functions for European style options are neither linear combinations of 
normal random variables nor they can be obtained by monotone transformations as for 
the case of the geometric average Asian options. To address the problem Imai and Tan 
propose to approximate an arbitrary function g, such that g~^ is the payoff function of 
a European derivative contract, with its first order Taylor expansion: 

.9(e) =.9(6)+^/ _Aez (36) 
^ oei 



e— e 



The approximated function is linear in the standard normal random vector Ae and 
we can rely on the same results obtained in the previous subsection. By considering 
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an arbitrary point of expansion, such as e = 0, we can derive the first column of the 
optimal orthogonal matrix A*. We can find the complete matrix by expanding g at 
different points and then run the optimization algorithm. 

Summarizing the optimization can be formulated as follow: 

max ) (37) 

A kGR" \dei e=ej 

subject to ||A.k|| = 1 and A*j • A.k = 0, j = 1, . . . , k — 1, k < n. 

Although equation (|32p provides an easy solution at each step, the correct procedure 
requires that A.k must be orthonormal to all the previous (and future) columns. This 
feature can be easily obtained by the Gram-Schmidt orthonormalization or even better 
by the QR method that is numerically stabler. As for the latter we must note that 
the QR method might return opposite matrices at different time steps (cosmetic sign 
adjustment). This does not affect the problem because the solution in equation ([32|) 
can be either with a positive and negative sign. Furthermore, we stress that it is not 
necessary to run the complete QR method at each step. Indeed, all the columns already 
calculated are orthogonal and we should use a "partial" QR method that considerably 
reduces the computational burden as it will be shown in the appendix. 

Imai and Tan set e\ — 0, 62 = (1, 0, . . . , 0), . . . , e'k = (1, 1, 1, . . . , 0, . . . , 0), . . . , = 
(1, , . . . , 1, 0)', the fc-th point has k-1 leading ones. 

The choice is arbitrary and a different set can be used that would return different 
optimal orthogonal matrices. 

Moreover the computational cost can be reduced by only seeking a suboptimal 
matrix with optimal columns up to k* < n. This approximation is reasonable since 
in practice only a few components are of relevance as will be shown in the numerical 
examples. 

4.3 Asian Options Case 

We consider the function g — g — K in equation (|14p . it is then, easy to verify that its 
variance can be expressed as: 



MN 



MN MN / 

{9 (e)) = ^ ^ ea:p + fi, + (1/2) ^ (Cf, + C 
1=1 j=i \ 



1=1 



exp 



(MN \ 



Due to the tractability of the function above, Imai and Tan provide some imple- 
mentations of the LT construction. We only show two of them. 

The variance contribution for the first p dimensions can be defined as: 



MN MN / p 

{9 (e)) = 5] 5] exp + /i, + (1/2) J2 {Cfi + C\ 
i=i j=i \ 



1=1 



exp 



(39) 

Working with algebra and approximating the exponential in the square bracket up to 
the first order, we can obtain the first formulation for the optimal matrix A: 

MNMN / P \ 

A ^^'L. = E E + + (1/2) E ((^^)' + (^^ f) CaC.i (40) 



i=l 3 = 1 



1=1 



subject to ||A.p|| = 1 and A.j • A*p — for j < k. 

The second formulation consists in applying the general approach by expanding the 
function of equation (|14p up to the first order: 



NM /NM I NM \ \ 

7(e) = g(e) + E E ^"^"P + E '-'^^^^ Ae/ 
1=1 \i=i \ k=i J J 



(41) 
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We start the optimization procedure by finding the first column of the optimal matrix 
A: 



NM /NM \ 

g{t) = g{0) + 51 H exp Ca Aei 

1=1 / 



(42) 



we set a,- = 



J2^Ji exp in,) Cu'j = Em=i {j2^Ji exp Cf^/^^ A„u in order to for- 
mulate equation ^ as ([28]). Set d^^) = (e^S . . . ,e^^^")'^ and B^^) = (d(i))'^C^'' we 
know from the linear combination case that A*^ = ± ugli) || ■ 

The p-th optimal column can be found considering the p-th starting point of the 
Imai and Tan's strategy. This results in: 

NM /NM / p-1 \ \ 

5(e) = 5(ep) + E E U» + E (43) 

1=1 \i=l \ k=l / / 

where C*)., k < p have been already found at the p — 1 previous steps and A.p 
must be orthogonal to all the other columns. As for the first step we define d'P^ — 

(^exp (a^i + E£} Cl,}j ,...,exp (^mw + EVi C^j^k) ^) and B(p) = {d(p))'^C^>^ the 

solution is A*p = ± ||b(p)|| • 

Alternatively, the optimal A* can be equivalently obtained by calculating the eigen- 
vector corresponding to the largest eigenvalues of the following matrix: 

NM NM / p-1 \ 

E E ^^^^ + A^:. + E + C:^) Cp'-C^'^ = dm5(d(P))EM^d*a5(d(P)) 

i=l 3 = 1 \ k=l j 

(44) 

and after imposing the orthonormality condition by the QR method. 



5 Simulation Framework 

We consider the constant volatility case only, and run our simulation with different 
combinations of path-generation techniques and different random number generators. 

As far as path-generation methods are concerned we use the standard Cholesky, 
the PCA and the two LT decompositions for Asian options introduced in the previous 
subsections. In particular for the first two approaches we rely on the properties of the 
Kronecker product in order to compute the decomposition fast (see Dahl, Benth [2] and 
[3] and Sabino [l3| for further details). 

LT methods require the iterative calculations of orthogonal matrices. We attain the 
task implementing an ad hoc QR factorization, as described in Appendix, that does not 
require high computational cost. For the LT decomposition the total computational 
time is than the sum of the time to compute the Cholesky and the optimal orthogonal 
matrix A. 

The numerical test consists of three main steps: 

1. Random number generation by standard MC, LHS or RQMC. 

2. Path generation with Cholesky, PCA, and the two LT algorithm discussed above 
(LTl and LT2, respectively). 

3. MC estimation. 

As RQMC generator we use a Faure-Tezuka scrambled version of the 50-dimensional 
Sobol sequence satisfying Sobols property A (see Glasserman [4], Jackel [9] and Owen 
|12j for further details). We pad the remaining random components out with LHS. This 
strategy is intended to investigate the effective improvement of the LT methods when 
coupled with QMC. Indeed, it can be proven that the LHS gives good variance reduc- 
tions when the target function is sum of one-dimensional functions (see Glasserman [3] 
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SiiO) 


= 100 




K 


= 90, 100 and 110 




r 


= 4% 




T 


= 1 




o-i 


= 10% + ^40% for i = 1, . . 


.,10 


Pij 


= and 40% forii,j = l,. 


..,10 



Table 1: Input Parameters. 



and Owen [H]). On the other hand, the LT methods is conceived to capture the lower 
effective dimension in superposition sense for hnear combinations. As a consequence, 
we should already observe a high accuracy when running the simulation with LHS and 
LT. Our setting is thought to test how large is the improvement given by the LT fac- 
torization. We compute a suboptimal A up to dimension 50 in order to be coherent 
with the choice of the 50-dimcnsional Sobol sequence. 

Stratification introduces correlation among random drawings so that the hypothe- 
sis of the Central Limit Theorem are not satisfied and we cannot compute the RMSE 
straightforward. We rely on the batch methods that consists of repeating Nb simula- 
tions for B times (batches). 

6 Numerical Investigations 

We develop our simulation procedure in order to test the computational burden and the 
efficiency of the Linear Transform method. We compare its results with those obtained 
with standard techniques like Cholesky and PCA decompositions. Furthermore, we use 
several random number generators, in particular, we adopt a Faure-Tezuka scrambled 
version of the 50-dimensional Sobol' sequence satisfying the Sobol's property A. 

As a numerical example, we estimate the fair price of an Asian option on a basket 
of M = 10 imder lying assets with N — 250 sampled points. 

The chosen parameters are those in the original paper of Imai and Tan ^ and are 
shown in Tabled! 

The nominal dimension of the problem is M x = 2500 equal to the number of 
rows and columns of the global correlation matrix Sa/jv- 

We perform the path-generation by computing the Cholesky, the PCA and two 
versions of the LT decompositions of the global correlation matrix S Af jv ■ We label LT2 
for the general case and LTl for the method described for Asian options only. As far 
as the first two generations are concerned, we rely on the properties of the Kronecker 
product in order to reduce the computational burden as described in Dahl, Benth [2] 
and [3^ and Sabino [T5] . 

As far as the implementation of the two LT methods proposed by Imai and Tan is 
concerned, we apply the fast version of the QR decomposition described in the appendix. 

The simulation procedure is implemented in MATLAB running on a laptop with an 
Intel Pentium M, processor 1.60 GHz and 1 GB RAM. 

Table [2] shows the percentage of the cumulative contribution of the variance for the 
first 10 components both for the zero and positive correlation cases, this is the ratio 
between equation ((55)) and ([55)) with p up to 10. All results up to p = 5 are consistent 
with those presented by Imai and Tan [5]. It can be noticed that the LT is the best 
performing path-generation technique in the statistical sense specified above, where the 
first specification is a bit better. The PCA decomposition is almost as accurate as the 
LT approach for the correlation case only. 

The effective dimensions found with each method are reported in Table [5| The 
advantage of the PCA and LT methods with respect to the Cholesky decomposition is 
evident both for the correlation and uncorrelation cases. The Cholesky decomposition 
collects 98.58% and 98.70% of the total variance for p = 2000 for the uncorrelation and 
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Uncorrelation 


Correlation 


Dimension 




PCA 


LTl 


LT2 




PCA 


LTl 


LT2 


1 




23 1 .5 


88 41 


88 41 


41 


91 40 


94.41 


94.41 


2 




42.24 


QO fi8 


90 24 


60 


93 04 


95 27 


95 1 7 


3 


n 07 




Q3 53 


92.20 


71 


94 32 


96 25 


95 91 


4 




fiQ Q1 


Q5 1 


94.11 


79 


95 28 


97 02 


96 58 

tyW. tJO 


5 


0.21 


7Q 30 


Qfi fiQ 


95 45 


85 


95 99 


97 59 


97 25 


6 


31 


8fi 25 


Q7 83 


9fi 68 


91 


97 42 


97 93 


97 58 

ty ( . <JO 


7 


45 


91.14 


Q8 32 


97 53 


95 


97 93 


98 34 


98 04 


8 


0.61 


94.33 


98.60 


98.10 


1.00 


98.28 


98.61 


98.24 


9 


0.81 


94.87 


98.85 


98.45 


1.04 


98.51 


98.75 


98.45 


10 


1.05 


95.28 


98.97 


98.64 


1.08 


98.78 


98.92 


98.61 



Table 2: Percentage of Variance Contribution up to dimension 10. 



Uncorrelation 


Correlation 


Cholesky PCA LTl LT2 


Cholesky PCA LTl LT2 


dr > 2000 dr = 18 dr = 11 dr = 13 


dr > 2000 dr = 12 dr = 12 dr = 12 



Table 3: Effective Dimensions. 



correlation cases, respectively. 

We compute the computational times elapsed to decompose the global covariance 
matrix with each method so that wc can compare the efhciency of all the methods; we 
compute only 50 optimal columns for the LT technique. Table |4] shows the estimated 
times in seconds. 

The computational times we found are a lot lower than those presented by Imai and 
Tan [5j despite the fact that are computed with a slower computer. In particular, the 
implementation of the LT method with the QR approach presented in the appendix (up 
to 50 columns) is more efficient of a factor tirthy. Furthermore, the LT methods has the 
versatility to allow the computation of a suboptimal matrix that is statistically justified 
by ANOVA considerations. In contrast, the PCA decomposition lacks this possibility 
without losing information. 

In the case of time-depending volatilities we could not rely on the properties of 
the Kronecker products in order to reduce the computational costs to run the PCA 
decomposition of the global covariance matrix. In contrast, the ad hoc QR approach 
for the LT method is preserved and needs a computational time of the same order as 
we will present in future studies. 

In the case of time-dependent volatilities it is fundamental to implement a fast 
Cholesky decomposition to be coupled with the QR method (see Sabino pj for further 
details of this type of Cholesky algorithm) . 

As a final step, we launch a MC simulation in order to estimate the fair price of the 
Asian basket option with 8192 generations and 10 replications. 

As already mentioned, we use a standard pseudo-random generator, the LHS method 
and a 50-dimensional Faure-Tezuka scrambled version of the Sobol' sequence satisfying 
the Sobol's property A. 





Uncorrelation 


Correlation 




Cholesky 


PCA 


LTl 


LT2 


Cholesky 


PCA 


LTl 


LT2 


time 


0.60 


25.77 


53.14 


53.21 


0.59 


25.55 


53.02 


53.20 



Table 4: Computational Times. 
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Standard MC 






K=90 


K=100 


K=110 








1 H.^c 


RMSE 




RMSE 


Cholesky 


12.257 


0.038 


5.604 


0.029 


2.007 


0.018 


PCA 


12.291 


0.038 


5.648 


0.029 


2.040 


0.018 


LT1=50 


12.240 


0.038 


5.681 


0.029 


2.004 


0.018 


LT2=50 


12.256 


0.038 


5.687 


0.029 


2.007 


0.018 






LHS 






Price 


RMSE 


Price 


RMSE 


Price 


RMSE 


Cholesky 


12.3320 


0.0097 


5.670 


0.013 


2.0394 


0.0084 


PCA 


12.3292 


0.0025 


5.6655 


0.0032 


2.0389 


0.0034 


LT1=50 


12.3291 


0.0015 


5.5968 


0.0019 


2.0332 


0.0022 


LT2=50 


12.2468 


0.0022 


5.6015 


0.0017 


2.0348 


0.0017 






RQMC 






Price 


RMSE 


Price 


RMSE 


Price 


RMSE 


Cholesky 


12.3410 


0.0094 


5.631 


0.014 


2.021 


0.012 


PCA 


12.32900 


0.00060 


5.65770 


0.00039 


2.03360 


0.00041 


LT1=50 


12.32800 


0.00036 


5.65720 


0.00040 


2.03400 


0.00021 


LT2=50 


12.32800 


0.00025 


5.65690 


0.00019 


2.03420 


0.00039 



Table 5: Correlation Case: Estimated Prices and Errors. 



It is known that (R)QMC simulations do not yield any improvements with respect 
to standard MC ones when the problem dimension is high (generally d > 20/30). Owen 
[TT] proposes mainly two approaches to extend the better convergence of the (R)QMC 
in high dimensions: the Latin Super Cube method and the padding with LHS. 

Briefly, the former consists of grouping the input variables and rearranging their 
order with a random permutation. The latter consists in fixing the more important 
variables and then pad the remaining ones out with the LHS. 

Even if this last method requires more computational costs, it can give further 
insight into the LT method. Indeed, it can test if the LT really selects the best variables 
in statistical sense and reduces the effective dimension. For the presented case we 
compare its results with those obtained with the pure LHS generator. We then, choose 
a 50-dimensional SoboF sequence, coherent with the suboptimal matrix A, and pad the 
remaining 2450 dimensions out. 

Tables [5] and m show the results of our numerical experiment. All values are statis- 
tically consistent, but exhibit a different accuracy. 

As expected, standard Cholesky decomposition is almost not sensitive to the used 
random generation technique and gives the worst results. 

LT and PCA decompositions provide good improvements for the RMSEs for both 
the RQMC and LHS generations. As far as the last method is concerned, we note 
that it is sensitive to the decomposition used and returns lower RMSEs when the LT 
decomposition is applied. This means that the LT approach is really reducing the 
effective dimension in superposition sense, "splitting" the integrand function into a 
sum of linear functions. 

As already mentioned, the LHS should reduce the RMSE in the case the integrand 
function is the sum of one-dimensional functions. This is best accomplished by the LT 
as evident from the above results. 

The RQMC simulation and the LT decompositions confirm their superior perfor- 
mance. 

It can be noted that the RQMC is sensitive to the used decomposition approach and 
does not have any advantage over the LHS when we use the Cholesky decomposition. 
Our evaluations return RMSEs with the same accuracy as Imai and Tan |^ when 
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Standard MC 






K=90 


K=100 


K=110 






± ^. 1 V 1 1 1 J 


PncG 


RMSE 




RMSE 


Cholesky 


11.560 


0.021 


3.426 


0.015 


3625 


0.0051 


PCA 


11.553 


0.021 


3.414 


0.015 


0.3591 


0.0052 


LT1=50 


11.454 


0.021 


3.356 


0.015 


0.3525 


0.0051 


LT2=50 


11.454 


0.021 


3.357 


0.015 


0.3523 


0.0051 






LHS 






Price 


RMSE 


Price 


RMSE 


Price 


RMSE 


Cholesky 


11.5915 


0.0037 


3.432 


0.007 


0.3605 


0.0038 


PCA 


11.5913 


0.0064 


3.4546 


0.0054 


0.3686 


0.0030 


LT1=50 


11.4754 


0.0013 


3.3666 


0.0023 


0.3662 


0.0021 


LT2=50 


11.4779 


0.0030 


3.3693 


0.0022 


0.3662 


0.0015 






RQMC 






Price 


RMSE 


Price 


RMSE 


Price 


RMSE 


Cholesky 


11.5890 


0.0033 


3.4351 


0.0049 


0.3605 


0.0035 


PCA 


11.59000 


0.00039 


3.4444 


0.0015 


0.3662 


0.0011 


LT1=50 


11.59100 


0.00025 


3.44360 


0.00039 


0.36673 


0.00034 


LT2=50 


11.59200 


0.00036 


3.44440 


0.00033 


0.36599 


0.00034 



Table 6: Uncorrelation Case: Estimated Prices and Errors. 



we only consider a 50-dimensional Sobol' sequence without using the complete LSS. 

Our framework is more extreme and the LT provides the same efficiency for all the 
strike prices and all correlations considered. In contrast, the PCA approach gives high 
improvements only in the correlation case. 

The general and the Asian options settings of the LT decompositions are almost 
equally performing with the latter one giving slightly better results. 

We can conclude that the LT is the best decomposition method and tremendously 
enhances QMC simulations because it optimally reduces the effective dimension of the 
problem. 

The LT construction can be made faster from the computational point of view, 
provided we implement the QR decomposition described in the appendix. 

7 Conclusion 

In this paper we investigate the accuracy of the LT, introduced by Imai and Tan, both 
from the computational and the accuracy points of view. In particular, we implement a 
numerical procedure based on the QR factorization that realize the LT decomposition 
fast. Moreover, we extensively investigate the improvements the LT gives to QMC 
methods that is sensitive to the effective dimension of the problem. 

As a numerical test we launch a high-dimensional simulation with the same set of 
parameters as in Imai and Tan [8] in order to price Asian basket options. 

Our setting is more extreme than the one discussed in the cited references. We do 
not rely on the complete LSS high-dimensional extension of the features of the QMC 
but we use a lower dimensional scrambled Sobol' sequence only, and pad the remaining 
ones out with LHS. 

We compare these results with those published by Imai and Tan [5] and those we 
found when using different decompositions and different random number generators. 

The LT construction provides the best accuracy with respect to the standard Cholesky 
approach and the PCA decomposition. 

It provides considerable improvements even when simulations are carried out with a 
partial RQMC method. The LT accuracy is still notably better than the one we found 
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with the complete LHS. In particular, we attain RMSEs of the same order as those 
presented by Imai and Tan. 

Moreover, the fast QR decomposition we implement gives an improvement of a 
factor 30 in terms of computational time compared to the results presented in Imai and 
Tan [5] calculated with a slower computer. 

PCA decomposition enhances QMC simulations but still requires a high computa- 
tional burden when time-dependent volatilities are considered and does not give the 
versatility to find a suboptimal matrix without introducing bias (see Sabino I13j for 
details) . 

Our QR-implementation makes the LT more efficient and computationally more 
convenient while maintaining its versatility for different problems. 
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8 Appendix 



8.1 The QR Method 

The QR factorization of an m-by-n matrix A is given by: 

A = QR 



(45) 



where Q € R™^™ is orthogonal and R e ]]j™xn jg upper triangular. A fundamental 
result is that if A has full column rank, then the first n columns of Q form an or- 
thonormal basis of ran(A). As a consequence, the QR factorization provides a way to 
return an orthonormal basis for a set of (independent) vectors. Different approaches 
can be chosen to calculate the QR decomposition such as the Householder and Givens 
transformations (sec Golub, Van Loan f5| as a fundamental reference). 

The former transformations are rank-two corrections of the identity of the form: 



1 ■ ■ 








■ ■ 
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• • 


COS 9 


sin 6 


.. 










.. 


• • 


— sin 6 ■ 


cos 9 . 


.. 










.. 


• • 








.. 1 




i 


k 





(46) 



(47) 

G{i, k, 9) performs a counterclockwise rotation of 9 radians in the (i, k) plane. 

Householder and Givens transformations are orthogonal transformations constructed 
in order to introduce zeros in a vector. Indeed, suppose we are given with 7^ x g M" 
we can find H and G, the Householder and Given Rotation, respectively, that annihilate 
the fc-th component of x: 



(48) 

y, yfc = (49) 

The following scheme illustrates the idea for QR factorization with Givens rotations: 



H = y,yk 
G^x 



A = 



XXX 

X X 
X X 
X 



(2,3)GfA = 



XXX 

X X 
X 
X 



(?,A)gIgIa^ R 



Here we have highlighted the 2- vectors that define the underlying Given rotations. 
Generally Q — Gi ■ ■ ■ G„ where n is the total number of rotations and R = Q^A. 

Consider A e R"^", its QR decomposition A = QaRa and B € R™x("+i) with 
the first n columns equal to A, Qb and Rb- The QR decomposition of B can be easily 
obtained from Qa and Ra- 

Denote b the last column of B, so that _B = [A b] it leads to Q^B = [Ra Q'a^] — 
Rb- Rb has the following form: 



Rb — 



X X X X 

X X X 

X X 

X X 

X X 



(50) 



In order to obtain the complete QR factorization of B wc only need to find t Givens 
transformations Gi,...,Gf that introduce zeros in the n-th column making Rb = 
GJ, . . . , GJRb upper triangular. 

Summarizing Qb = Gi, . . . , GiQa and Rb ^ Gj, . . . , GJQaB. 
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